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A  two-phase  non-isothermal  model  is  developed  to  explore  the  interaction  between  heat  and  water 
transport  phenomena  in  a  PEM  fuel  cell.  The  numerical  model  is  a  two-dimensional  simulation  of  the 
two-phase  flow  using  multiphase  mixture  formulation  in  a  single-domain  approach.  For  this  purpose,  a 
comparison  between  non-isothermal  and  isothermal  fuel  cell  models  for  inlet  oxidant  streams  at  different 
humidity  levels  is  made.  Numerical  results  reveal  that  the  temperature  distribution  would  affect  the 
water  transport  through  liquid  saturation  amount  generated  and  its  location,  where  at  the  voltage  of 
0.55  V,  the  maximum  temperature  difference  is  3.7  °C.  At  low  relative  humidity  of  cathode,  the  average 
liquid  saturation  is  higher  and  the  liquid  free  space  is  smaller  for  the  isothermal  compared  with  the  non- 
isothermal  model.  When  the  inlet  cathode  is  fully  humidified,  the  phase  change  will  appear  at  the  full  face 
of  cathode  GDL  layer,  whereas  the  maximum  liquid  saturation  is  higher  for  the  isothermal  model.  Also, 
heat  release  due  to  condensation  of  water  vapor  and  vapor-phase  diffusion  which  provide  a  mechanism 
for  heat  removal  from  the  cell,  affect  the  temperature  distribution.  Instead  in  the  two-phase  zone,  water 
transport  via  vapor-phase  diffusion  due  to  the  temperature  gradient.  The  results  are  in  good  agreement 
with  the  previous  theoretical  works  done,  and  validated  by  the  available  experimental  data. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Water  management  is  extremely  important  for  balancing  the 
operation  of  PEM  fuel  cells,  to  avoid  flooding  while  maintaining 
proper  membrane  hydration  to  achieve  the  best  possible  perfor¬ 
mance.  The  membrane  in  PEM  fuel  cell  must  be  fully  hydrated  to 
get  the  optimal  proton  conductivity.  However,  due  to  low  operat¬ 
ing  temperatures  (70-90  °C),  PEM  fuel  cells  are  prone  to  gas-liquid 
formation,  particularly  when  they  are  highly  humidified  or  at  low 
gas  flow  rate  conditions  [1,2]. 

When  GDL  and  catalyst  layer  become  saturated  with  water 
vapor,  the  produced  water  vapor  starts  to  condense  and  therefore 
blocks  the  open  pores,  reducing  the  oxygen  transport  to  catalyst 
layer.  Flooding  becomes  a  major  factor  limiting  the  PEM  fuel  cell 
performance.  Thermal  management  is  also  required  to  remove  the 
excessive  heat  produced  due  to  various  heat  generation,  including 
the  irreversible  heat  from  electrochemical  reactions,  entropic  heat, 
Joule  heating  arising  from  the  electrolyte  ionic  resistance  and  the 
heat  from  condensation  that  could  dry  out  the  membrane.  These 
heat  sources,  especially  the  irreversible  reaction  heat  and  entropic 
heat  could  rise  the  fuel  cell  temperature  during  operation  and  insuf¬ 
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ficient  cooling  may  result  in  local  hot  spots  which  could  reduce 
hydration,  therefore  hampering  membrane  performance  [3].  Also 
having  a  uniform  temperature  distribution  in  the  porous  electrode 
and  small  temperature  variation  is  favored  where  proton  conduc¬ 
tivity  of  the  polymer  electrolyte  membrane  strongly  depends  on 
the  degree  of  its  hydration  which  is  affected  by  temperature.  In  this 
respect,  a  numerical  model  is  a  useful  tool  for  understanding  water 
and  thermal  management  and  interaction  between  them  that  is 
essential  for  proper  PEM  fuel  cell  operation. 

A  number  of  different  computational  approaches  for  PEM  fuel 
cells  have  been  carried  out  in  recent  years  [1,3-23],  considering 
two-phase  or  single-phase  water  transport  with  or  without  heat 
transfer  at  proton  exchange  membrane  regions.  Several  studies 
[4-7]  have  attempted  to  predict  the  temperature  distribution  for 
single-phase  condition.  Ju  et  al.  [8]  reviewed  the  single-phase,  non- 
isothermal  models  of  PEM  fuel  cells  in  detail.  Two-phase  transport 
in  PEM  fuel  cells  has  also  been  studied  by  several  researchers 
[1,9-12];  however,  the  focus  of  these  studies  was  primarily  on 
the  isothermal  investigation  of  the  transport  phenomena.  The 
two-phase  non-isothermal  model,  gives  a  proper  simultaneous 
description  of  water  and  thermal  management  with  phase  change. 

Rowe  and  Li  [13]  and  Mishra  et  al.  [14]  performed  a  study  on 
water  and  thermal  management  in  PEM  fuel  cell  considering  a 
steady-state,  one-dimensional  approach.  Berning  and  Djlali  [15] 
presented  a  model  based  on  the  unsaturated  flow  theory  with 
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Nomenclature 

a  water  activity  or  specific  electrochemically  active 

area  (m2  m-3) 

A  superficial  electrode  area  (m2) 

C  molar  concentration  of  species  i  (mol  m-3 ) 

cp  specific  heat  capacity  (J  kg-1  K-1 ) 

Dl  mass  diffusivity  of  species  i  (m2  s-1 ) 

EW  equivalent  weight  of  dry  membrane  (kg  mol-1 ) 

F  Faraday  constant,  96,487  (C  mol-1 ) 

hfg  latent  heat  of  evaporation  of  water  (kj  kg-1 ) 

i'o  exchange  current  density  (A  m-2 ) 

I  current  density  (A  m-2 ) 

j  transfer  current  (A  m-3 ) 

j  Leverett  function 

j\  mass  flux  of  liquid  phase  (kg  m-2  s-1 ) 

k  thermal  conductivity  (W  nrr 1  K-1 ) 

/<rk  relative  permeability  of  phase  k 

I<  hydraulic  permeability  (m2 ) 

M1  molecular  weight  of  species  i  (kg  mol-1 ) 
m/j'  mass  fraction  of  species  i  in  phase  k 
/hfg  mass  flow  rate  (kgs-1) 

n  the  direction  normal  to  the  surface 

nd  electro-osmotic  drag  coefficient 

P  pressure  (Pa) 

R  universal  gas  constant,  8.314  J  mol-1  K-1 

RH  relative  humidity 

s  liquid  saturation 

S  source  term 

T  temperature  (K) 

u  velocity  vector  (ms-1) 

Uq  thermodynamic  equilibrium  potential  (V) 

Vce\\  cell  potential  (V) 

V  volume  (m3) 

x  mole  fraction 

Greek  symbols 

a  transfer  coefficient  for  reaction 

s  volume  fraction 

£mc  volume  fraction  of  ionomer  phase  in  catalyst  layer 
0  potential  (V) 

yc  advection  coefficient  for  transport  of  species 
yh  advection  coefficient  for  heat  transfer 

I I  overpotential  (V) 

k  ionic  conductivity  (S  m-1 ) 

X  polymer  water  content 

Xk  relative  mobility  of  phase  k 

fi  viscosity  (Pas) 

v  kinematic  viscosity  (m2  s-1 ) 

0c  contact  angle  (°) 

p  density  (kg  m-3) 

Pdry.m  dry  membrane  density  (kg  m-3 ) 
a  surface  tension  (Nm-1 ) 

£  stoichiometry  flow  ratio 

Subscripts 
a  anode 

c  cathode  or  capillary 

e  electrolyte 

in  inlet 

g  gas 

1  liquid 

ref  reference 

sat  saturation 


s  solid 

0  standard  condition,  298.15  K  and  101.3  kPa 

Superscripts 
eff  effective 

H2  hydrogen 

H20  water 

02  oxygen 


the  assumption  of  constant  gas  pressure  across  the  porous  media. 
In  this  model  two  separate  computational  domains  had  to  be  set 
up,  one  for  gas  flow  channels  and  the  other  for  GDL  to  accom¬ 
modate  the  heat  transfer  through  the  solid  matrix  of  the  porous 
media.  Mazumder  and  Cole  [16]  adopted  the  M2  model  to  inves¬ 
tigate  the  formation  and  transport  of  liquid  water  in  PEM  fuel 
cells  for  a  three-dimensional  geometry.  In  this  model  the  cell  was 
assumed  to  be  isothermal  by  setting  the  thermal  conductivity  in  all 
regions  to  be  very  high,  therefore  neglecting  the  thermal  effects  on 
phase  change.  Birgersson  et  al.  [17]  presented  a  two-dimensional 
model  based  on  the  multifluid  approach,  which  is  in  contrast  to 
the  M2  model.  Wang  and  Wang  [18]  developed  a  non-isothermal, 
two-phase  model  based  on  the  M2  approach  and  identified  the 
importance  of  water  transport  as  well  as  heat  removal  via  vapor- 
phase  diffusion  with  variable  temperature.  The  model  presented 
by  Hwang  [19]  illustrates  the  behaviors  of  the  two-phase  flow  and 
heat  transfer  in  a  porous  electrode.  Pasaogullari  et  al.  [20]  devel¬ 
oped  a  model  to  investigate  heat  and  mass  transfer  in  the  GDL 
cathode  simultaneously.  Also  Yuan  and  Sunden  [21]  presented  a 
non-isothermal,  two-phase  model  in  multi-dimensional  situations. 
However,  only  the  cathode  GDL  and  gas  channel  were  considered. 
Afshari  and  Jazayeri  [22,23]  developed  a  mathematical  model  for 
the  heat  transfer  and  liquid  water  formation  in  a  PEM  fuel  cell 
and  investigated  the  thermal  and  water  management  effects  on  cell 
performance. 

In  the  present  study,  a  two-dimensional,  two-phase,  non- 
isothermal,  and  single-domain  model  together  with  coupled 
electrochemical  relations  is  analyzed  and  the  effects  of  interactions 
between  water  and  heat  transfer  with  phase  changes  in  a  PEM  fuel 
cell  are  investigated. 


2.  Mathematical  modeling 


The  model  domain  consists  of  the  following  subregions:  the  gas 
channels,  gas  diffusion  layers  (GDLs),  catalyst  layers  for  both  the 
anode  and  cathode  sides,  with  the  membrane  in  the  middle.  Fuel 
and  oxidant  flow  through  channels  and  are  distributed  into  anode 
and  cathode.  Hydrogen  oxidation  and  oxygen  reduction  reactions 
occur  only  within  the  active  catalyst  layers.  The  fuel  and  oxidant 
flow  rates  can  be  described  by  a  stoichiometric  flow  ratio,  £ ,  defined 
as  the  amount  of  reactant  in  the  gas  chamber  fed  to  the  amount 
required  by  the  electrochemical  reaction  [18].  That  is 


— 


0-l2,in^a,in 


2F  Aa  jn 
/ref  A 


(1) 


= 


Q)2,in^c,ir 


4F  Ac  jn 
/ref  A 


(2) 


where  A  is  the  superficial  electrode  area,  and  Aa,  Ac,  Din  a  and  U[n  c 
are  the  flow  cross-sectional  areas  and  the  inlet  velocities  of  the 
anode  and  cathode  gas  channels,  respectively.  The  stoichiometric 
flow  ratios  defined  in  the  present  work  for  Eqs.  (1 )  and  (2)  are  cho¬ 
sen  at  the  fixed  reference  current  density  of  1 A  cm-2.  Therefore  the 
flow  rates  of  fuel  and  oxidant  are  constant. 
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2.1  Assumptions 

This  model  assumes: 

•  Two-dimensional. 

•  Steady-state  flow. 

•  Ideal  gas  mixtures. 

•  Laminar  and  incompressible  flow  with  low  Reynolds  numbers 
and  pressure  gradients. 

•  Negligible  ohmic  drop  in  the  electronically  conductive  solid 
matrix  of  porous  electrodes,  and  catalyst  layers. 

•  Negligible  Soret,  Dufour,  gravity  and  radiation  effects. 

•  Isotropic  and  homogeneous  electrodes,  catalyst  layers  and  mem¬ 
brane  are  characterized  by  effective  porosities  and  permeabilities. 

•  Constant  viscosity  of  gas  mixture  is  calculated  for  the  inlet  con¬ 
dition. 

•  No  contact  resistance  at  the  interfaces  between  different  layers. 
2.2.  Governing  equations 

In  contrast  to  the  usual  approach  where  separate  differential 
equations  are  employed  for  different  regions,  here  a  unified  single¬ 
domain  approach  with  a  single  set  of  governing  equations  is  applied 
to  all  regions.  Based  on  the  multi-domain  method,  the  compu¬ 
tational  domain  is  divided  into  a  number  of  sub-domains  and 
different  sets  of  conservation  equations  are  developed  for  differ¬ 
ent  sub-domains  and  the  interfacial  boundary  conditions  establish 
the  connection  between  these  equations.  In  the  developed  single¬ 
domain  method,  one  set  of  conservation  equations  is  considered 
for  different  regions  of  a  PEM  fuel  cell.  As  a  result,  special  numeri¬ 
cal  treatments  have  been  used  by  defining  extremely  large  or  small 
physical  and  transport  parameters  in  a  region  [1,11  ].  Therefore  this 
model  has  a  set  of  coupled  non-linear  partial  differential  equa¬ 
tions  including  conservations  of  mass,  momentum,  species,  energy 
and  charge  with  electrochemical  relations.  The  multiphase  mix¬ 
ture  model  (M2)  is  used  for  the  two-phase  flow  that  is  an  exact 
formulation  of  classical  two-fluid,  two-phase  model  having  a  sin¬ 
gle  equation.  The  main  difference  between  the  M2  model  and  the 
unsaturated  flow  theory  is  that  the  former  does  not  require  the 
approximation  of  constant  gas-phase  pressure.  Another  important 
feature  of  the  M2  model  is  that  it  can  be  easily  used  in  a  prob¬ 
lem  where  single  and  double  phase  zones  co-exist  leading  to  a 
substantial  decrease  in  the  numerical  complexity  [1,24-29]. 

Mass  conservation:  Conservation  of  mass  for  the  two-phase  mix¬ 
ture  using  the  M2  model  is  as  follows: 


Momentum  conservation :  Conservation  of  momentum  for  the 
two-phase  mixture  with  mixture  velocity  u  can  be  written  as: 


-y  V  •( puu)  =  -Vp  + V  •  (/xVu)  +  SDa 


where  n  is  mixture  viscosity  defined  as  [24]: 

kr\  /<rg 
fl  =  P  —  +  — - 
V\  vg 


(5) 


(6) 


In  porous  regions,  superficial  velocities  are  used  in  order  to  auto¬ 
matically  ensure  mass  flux  continuity  at  the  interfaces  between 
porous  and  non-porous  regions.  Also  intrinsic  transport  properties 
in  the  porous  regions  are  changed  into  effective  transport  prop¬ 
erties  taking  into  account  the  effects  of  porosity  and  tortuosity 
using  Bruggeman  correlation  [1].  The  source  term  of  momentum 
equations  in  Table  1  is  employed  to  consider  Darcys  law  under  the 
limiting  condition  where  the  permeability  of  porous  media  is  small, 
resulting  in  low  velocity. 

Species  conservation:  Conservation  of  species  equation  using  the 
M2  model,  in  terms  of  molar  concentration  is  as  follows  [26]: 


Viri-uC)  =  V[Dj;effvq]  -  V 


C1 

Zl 

pg 


j\ 


+  5k 


(7) 


where  C  is  the  total  concentration  of  species  i  in  the  liquid  and 
gas  phases.  The  liquid  and  gas  phases  have  different  flow-fields  so 
the  advective  transport  of  species  is  corrected  using  the  following 
correction  factor,  yc  [1]: 


Yc  = 


P  ( 

CH2°  l  MH2° 

PK g 

PgO  -5) 


+  Ag 


Pg  ) 


for  water 
for  other  species 


(8) 


where  X\  and  Ag  are  the  relative  mobility  of  liquid  and  gas  phases, 
respectively  [24]. 


•\  _  ^rk/^k 

k'£iA,<M 


O) 


vk  is  the  kinematic  viscosity,  k  denotes  the  liquid  or  gas  phase, 
kT i  and  kTg  are  the  relative  permeabilities  of  liquid  and  gas  phases, 
respectively  [15]. 

krg  =  (l-S)3,  krl=S3  (10) 


V  ■  (pu)  =  0  (3) 

where  u  and  p  are  mixture  velocity  and  density,  respectively  [24,25  ]. 

p  =  pxs  + pg{\-s)  (4) 

s  and  (1  -s)  are  the  fraction  of  open  pore  space  occupied  by  the 
liquid  and  gas  phases,  respectively. 


The  pmf  is  mass  fraction  of  species  i,  in  terms  of  s: 

pmf  =  pgmPg{\  -  s)  +  pxmf{s  (11) 

The  second  term  on  the  right-hand  side  of  the  species  conserva¬ 
tion  equation  represents  the  capillary  transport  using  the  theory  of 
capillary  transport  in  hydrophobic  GDL  developed  by  Pasaogullari 


Table  1 

Source  terms  for  momentum,  species,  energy  and  charge  conservation  equations  for  various  regions. 


Equation 

Source  terms 

Flow  channels 

GDLs 

Catalyst  layers 

Membrane 

Momentum 

•Soar  =  0 

Soar  =  —  y  ^ 

Soar  =  —  y  ^ 

Soar  =  —  y  ^ 

Species 

sk=0 

sk  =  0 

on 

II 

1 

< 

1? 

1 

Sk  =  -V.  (£/) 

Energy 

ST  =  /ifgrhfg 

ST  =  hfgriifg 

St  -J  (jl  +  )  +  jiff  +  hfg^fg 

ST  =  +  hfgriifg 

Charge 

Se  =  0 

Se=  0 

Se=j 

Se=0 

(Mk: 

chemical  formula  of  species  k 

Electrochemical  reaction:  ^skMk  =  ne~  where  <  sk  :  stoichiometry  coefficient 

^  n  :  number  of  electrons  transferred 

anode  hydrogen  oxidation  reaction  (HOR)  H2  -  2H+  ->  2e_,  cathode  oxygen  reduction  reaction  (ORR)  2H20  -  02  -  4H+  ->  4e_. 
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Table  2 

Physical  and  transport  properties. 


Parameter 


Value 


Water  vapor  activity 

Ionic  conductivity  of  membrane  (S  m-1 ) 

Polymer  water  content 

Electro-osmotic  drag  coefficient 

Water  diffusivity  in  membrane  (m2  s_1 ) 

Gas  diffusivity  (m2  s-1 ) 

Gas  density  (kg  nrr3) 


“  ~~  Aat  “  cH20 

k  =  (0.5139A  -  0.326)exp  (l268.0  (  ^  -  I )  ) 

'  0.043  +  17.81a  -  39.85a2  +  36.0a3  0  <  a  <  1 
X  =  4  14  + 1.4(a  -  1)  l<a<3 

16.8  a  >  3 

1.0  X  <14 


n<d  ’  -^(A  -  14)  +  1.0  otherwise 

h2o  _  /  3.1  x  10“7A(eom  -  1)  e(-2346/T)  o  <  X  <  3 
e  “  \  4.17  xlO“8A(  1  +  161  e~x)  e^2346/T)  A  >3 

pg  =  — 

RT  \  (Wj/M;) 


and  Wang  [1],  In  the  absence  of  gravity,  liquid  flux,  ji  is  given  by 
[24,26]: 


Z+g 


KVpc 


(12) 


pc  is  the  capillary  pressure  defined  as  [25]: 

/  e  \  i/2 

pc  =  a  cos(0c)(-J  J(s )  (13) 

where  J(s),  Leverett  function  is  given  for  both  the  hydrophobic  and 
hydrophilic  GDL  layers  [30], 


J(s)= 


1.417(1  -  s)- 2.120(1  -s)2  +  1.263(1  -s)3,  if+  <  90° 
1.417s  -  2.120s2  +  1.263s3  if0c  >  90° 

(14) 


The  gas  phase  diffusion  coefficient  is  modified  to  be  effective 
using  Bruggeman  correlation  [  18  ]  to  account  for  the  effects  of  poros¬ 
ity  and  tortuosity  of  the  porous  electrodes  and  catalyst  layers.  Also 
in  the  two-phase  zones  where  the  gas  phase  could  be  saturated 
with  water  vapor,  the  equation  is  modified  by  the  liquid  saturation. 

Note  that  in  an  isothermal  model  the  water  vapor  concentra¬ 
tion  in  the  two-phase  zone  takes  the  saturation  value,  thus  the 
molecular  diffusion  of  water  vapor  disappear  (Dgeff  =  0)  [1  ].  In  the 
non-isothermal  situation,  the  water  vapor  saturation  concentration 
that  is  a  strong  function  of  temperature,  varies  and  therefore  the 
vapor  phase  diffusion  emerges  as  a  new  transport  mechanism  in 
the  two-phase  zone  [18]. 

+eff  =  [e(l  -s+'+j,  (15) 

The  sources  terms  of  species  conservation  equations  in  Table  1, 
are  the  representative  of  source  or  sink  terms  due  to  the  electro- 
osmotic  drag  of  water  in  the  membrane  and  electrochemical 
reactions  in  catalyst  layers. 

Electronic  charge  conservation :  Conservation  of  electronic  charge 
equation  is: 

V.(KeffV(/>e)  +  Se  =  0  (16) 


This  equation  describes  proton  transport  inside  the  mem¬ 
brane-electrode  assembly.  The  dependence  of  proton  conductiv¬ 
ity  on  water  content  is  calculated  using  the  empirical  expression 
of  Springer  et  al.  [31,8]  as  shown  in  Table  2.  The  effective  proton 
conductivity  in  porous  media  is  described  using  the  Bruggeman 
relation  [8]. 

/ceff  =  £1,5/c  (17) 

The  source  term  for  the  charge  equation  in  Table  1  represents 
the  transfer  of  current  between  the  electronically  conductive  solid 
matrix  and  the  electrolyte  for  both  the  anode  and  cathode  catalyst 
layers.  In  a  PEM  fuel  cell,  the  anode  hydrogen  oxidation  reaction 
exhibits  fast  electrokinetics  with  a  low  surface  overpotential;  there¬ 
fore,  it  can  easily  be  expressed  by  a  linear  equation,  whereas  the 
cathode  oxygen  reduction  reaction  has  a  relatively  slow  kinetics 
with  a  higher  surface  overpotential,  which  is  adequately  described 
by  Tafel  kinetics  and  is  summarized  in  Table  3. 

Energy  conservation :  Conservation  of  energy  equation  using  the 
M2  model  is  as  follow  [18]: 

W.(YhPCpUT)  =  V.(keffVT)  +  ST  (18) 

The  heat  terms  in  the  above-mentioned  equation,  as  shown  in 
Table  1,  contain  irreversible  heat  of  the  electrochemical  reaction, 
reversible  entropic  heat,  Joule  heating,  and  an  extra  source  due  to 
condensation  and  evaporation.  In  the  last  source  term,  hfg  is  the 
latent  heat  of  condensation  or  evaporation  and  rhfg  is  the  mass  flow 
rate  due  to  phase  change  that  can  readily  be  calculated  from  the 
continuity  equation  of  the  liquid  phase  [18,32]: 

m^  =  V-(/olUi)  =  V-(/1  +  A.1/ou)  (19) 

Since  the  pore  size  and  advection  in  the  GDL  are  relatively  small, 
the  temperature  in  solid  and  fluid  phases  are  assumed  to  be  equal. 
The  advective  term  in  Eq.  (18)  due  to  different  flow  fields  of  liquid 
and  gas  phases  is  corrected  using  the  following  correction  factor, 


Table  3 

Electrochemical  properties. 


Parameter 

Anode 

Cathode 

Transfer  current  density  (A  rrr3 ) 

Surface  overpotential  (V) 

Open  circuit  potential  (V) 

Transfer  coefficient 

Exchange  current  density  times  reaction  surface  area  (A  rrr3 ) 

l-0,a(|)  7  ( 

T]  =  0S  -  0e  -  00  (0s 

u0=o 

aa  +  cue  =  2 
aio,e  =  10  x  109 

=  0) 

j  =  -fli0,c  (^)  exp  (-IfFtj) 

0  =0S  -0e  -  Oo(0s  =  Vcell) 

U0  =  1.23  -  0.9  X  10-3  (T -  298.15) 

Qfc  =  1 

ai0,e(353)=2.0  x  104 

ch'o,c(t)  =  (1  -s)ai'o,c(353)  exp  (-16,456  (j  -  sfc))- 
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anode  channel  cathode  channel 


X 


Fig.  1.  Two-dimensional  geometry  and  computational  mesh  for  PEM  fuel  cell  model. 


Yh  [18] 

_  p(^lcp,l  +  ^gCp.g) 
sAcp,i  +  (1  -  s)PgCP,g 


(20) 


A  constant  temperature  is  applied  to  the  anode  and  cathode  gas 
inlets. 

^a.in  —  ^cell*  ^c,in  —  ^cell  (24) 


The  physical  and  transport  properties  in  the  above  mentioned 
equations  are  presented  in  Table  2. 

2.3.  Boundary  conditions 

The  complete  set  of  governing  equations  representing  the  math¬ 
ematical  model  is  given  by  Eqs.  (3),  (5),  (7),  (16)  and  (18)  that  forms 
a  set  of  equations  with  seven  unknowns:  u,p,  CH2,  C°2,  CH2°,  0e 
and  T.  The  advantage  of  this  model  is  that  only  the  external  bound¬ 
ary  conditions  are  required  in  this  single-domain  formulation.  The 
boundary  conditions  are  specified  as  follows. 


2.3.1.  Inlet  boundaries 

The  inlet  velocity  in  the  gas  channels  are  expressed  by  the 
respective  stoichiometric  flow  ratio: 


(21) 


Uain  and  Uc[n  are  derived  from  Eqs.  (1)  and  (2).  The  molar  con¬ 
centrations  of  species  are  determined  by  the  inlet  pressure  and 
humidity  according  to  the  ideal  gas  law. 


CH2°  =  yH?  Q  ^a'in 
'“'a, in  Aa,in  p t  . 

1X1  a, in 


=  RH, 


Psat 

^a,in 


^a,in 

^Ta,in 


Ch2°  —  yH?Q  ^c,in 

L'c,in  Ac,in  p t  . 

c,in 


=  RH, 


Psat  Pc, in 
C'mP^PT~ 


(22) 


Since  summation  of  mole  fraction  species  at  anode  and  cathode 
inlets  are  unit,  also  the  ratio  of  nitrogen  and  oxygen  in  dry  air  is 
known  to  be  79:21,  the  inlet  hydrogen  and  oxygen  concentration 
can  be  found  via: 


CH-2  =  (1  —  XH2°) — — 
a, in  v  a, in  ' 


a, in 


1  a, in 


C°2  —  y°2  ^c’in 
c,in  c,in  dt  . 

JVJc,ir 


(23) 


2.3.2.  Outlet  boundaries 

The  flow  at  the  outlet  is  assumed  to  be  fully  developed  or  no-flux 
with  a  given  back  pressure. 


(25) 


2.3.3.  Walls 

No-slip  and  impermeable  velocity  condition  and  no-flux  condi¬ 
tion  (except  at  the  outer  boundary  of  anode  and  cathode  channels 
for  thermal  boundary  conditions)  are  applied: 


S  =  0;  f-o;  f-o:  ”.o 

dn  on  on  on 


(26) 


At  the  outer  boundary  of  anode  and  cathode  channels,  either 
constant  temperature  or  constant  heat  flux  boundary  conditions 
can  be  imposed  [33].  Here  constant  temperature  boundary  condi¬ 
tion  is  used  [18]. 


^channel  wall  —  ^cell 


(27) 


2.4.  Numerical  procedures 

The  geometry  and  grid  structure  is  given  in  Fig.  1  with  the  spec¬ 
ifications  listed  in  Table  4.  Although  commercial  CFD  programs  are 
favored  and  widely  used,  still  in-house  programming  has  its  own 
advantages  due  to  a  better  control  of  a  detailed  modeling.  The  gov¬ 
erning  equations  were  discretized  using  a  finite-volume  method 
and  solved  using  a  developed  code  written  in  Fortran  software.  In 
this  code  the  pressure  and  velocity  fields  are  treated  using  SIMPLE 
pressure  correction  algorithm  for  a  single-domain  model.  It  should 
be  mentioned  that  in  spite  of  the  absence  of  some  species  partic¬ 
ularly  in  certain  regions  of  fuel  cell,  the  species  transport  equation 
can  still  be  applied  throughout  the  entire  computational  domain 
[34]. 

Stringent  numerical  tests  were  carried  out  to  ensure  that  the 
solution  was  independent  of  grid  size.  At  least  250  computational 
volumes  in  the  along  channel  and  50, 30, 10  and  60  computational 
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Table  4 

Dimensional  parameters  and  transport  properties. 


Description 

Unit 

Value 

Dimensional  parameters 

Cell  length 

mm 

70 

Channel  thickness 

mm 

1 

Gas  diffusion  layer  thickness 

fxm 

300 

Catalyst  thickness 

fxm 

10 

Membrane  thickness 

fxm 

50.8 

Operating  conditions 

Cell  temperature 

K 

353.0 

Anode/cathode  pressure 

atm 

1.5 

Anode/cathode  stoichiometry,  £ 

- 

2.0 

Anode  dry  gas  mole  fraction 

- 

0.0 

Cathode  dry  gas  mole  fraction 

- 

3.76 

Transport  parameters 

Reference  hydrogen  molar  concentration,  CJj^ 

mol  m3 

40.88 

Reference  oxygen  molar  concentration,  CjJ 

mol  m3 

40.88 

Faraday’s  constant,  F 

Cmol-1 

96,487 

Universal  gas  constant,  R 

J  mol-1  K-1 

8.34 

Anode  gas  viscosity 

pas 

1.101  x  10“5 

Cathode  gas  viscosity 

pas 

1.881  x  10“5 

Liquid  water  viscosity  at  80  °C 

pas 

3.56  x  10-4 

H2  diffusivity  in  the  anode  gas  channel 

m2  s-1 

5.457  x  10“5 

H20  diffusivity  in  the  anode  gas  channel 

m2  s-1 

5.457  x  10“5 

02  diffusivity  in  the  cathode  gas  channel 

m2  s_1 

1.806  x  10-5 

H20  diffusivity  in  the  cathode  gas  channel 

m2  s_1 

2.236  x  10“5 

Anode/cathode  GDL  porosity,  £Gdl 

- 

0.6 

Anode/cathode  catalyst  layer  porosity,  eCat 

- 

0.6 

Volume  fraction  of  membrane  in  catalyst  layers,  £mc 

- 

0.26 

Anode/cathode  GDL  permeability,  JCG dl 

m2 

6.875  x  10-13 

Contact  angle  of  GDL,  0C 

o 

110 

Surface  tension,  cr 

N  m1 

0.0625 

Equivalent  weight  of  membrane,  EW 

kg  mol-1 

1.1 

Dry  density  of  membrane,  pdry 

kg  m-3 

1.98  x  103 

Thermal  conductivity  of  hydrogen 

Witt1  K1 

0.2040 

Thermal  conductivity  of  oxygen 

Witt1  K-1 

0.0296 

Thermal  conductivity  of  nitrogen 

Wirr1  K-1 

0.0293 

Thermal  conductivity  of  water  vapor 

Witt1  K1 

0.0237 

Thermal  conductivity  of  liquid  water 

Witt1  K”1 

0.67 

Thermal  conductivity  of  membrane 

Wm1  K-1 

0.950 

Thermal  conductivity  of  GDL 

Witt1  K1 

0.5 

volumes  along  the  thickness  of  each  channel,  GDL,  catalyst  layer 
and  membrane,  respectively,  are  used.  For  a  two-dimensional  flow 
field,  this  grid  requirement  results  in  approximately  60,000  com¬ 
putational  cells.  The  coupled  set  of  equations  was  solved  iteratively, 
and  the  solution  was  checked  to  be  convergent  when  the  relative 
error  in  each  field  between  two  consecutive  iterations  was  less  than 
10-6.  The  flowchart  for  this  simulation  is  shown  in  Fig.  2.  The  CPU 
time  range  is  between  5  and  6  h  on  a  Pentium  IV  PC  (3.2  GHz,  1  GB 
RAM). 

3.  Results  and  discussion 

3.1  Model  validation 

For  validation  purpose,  the  data  of  a  fuel  cell  operated  under 
no  humidity  and  low  operating  temperature  (50  °C)  is  considered 
[35].  The  simulation  is  carried  out  using  the  developed  model  for 
the  same  experimental  conditions.  The  derived  polarization  curve 
is  compared  with  both  the  experimental  data  presented  by  Buchi 
and  Srinivasan  and  the  single-phase  model  presented  by  Um  and 
Wang  [10]  (see  Fig.  3).  The  obtained  results  from  this  model  show 
very  good  agreement  with  the  measured  experimental  data  with  a 
deviation  of  less  than  7%  which  could  be  as  a  result  of  2D  modeling 
where  the  rib  effects  are  not  considered.  At  high  current  densi¬ 
ties,  the  results  are  in  better  agreement  with  the  experimental  data 
compared  with  Um  and  Wang  results.  This  is  due  to  the  fact  that 
their  model  is  a  single  phase  one  where  the  effects  of  phase  changes 


Fig.  2.  Flowchart  of  the  simulation  process. 


at  high  current  densities,  specially  phase  changes  at  low  operating 
temperature  are  neglected. 

3.2.  Effect  of  temperature  on  phase  change 

In  this  section,  a  two-phase  non-isothermal  model  is  developed 
to  simulate  the  two-phase  phenomena  in  a  single  straight-channel 
of  a  PEM  fuel  cell.  Special  attention  has  been  given  to  the  impact 
of  the  interfacial  liquid  accumulation  on  GDL,  temperature  distri¬ 
bution,  and  interaction  between  phase  change  and  heat  transfer. 


Fig.  3.  Comparison  presented  study  with  an  experimental  polarization  curve  and  a 
single-phase  model. 
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Fig.  4.  (a)  Liquid  saturation  distribution  at  cathode  GDL  for  non-isothermal  model 
at  Vceii  =  0.55  V  for  AFCL  case  (b).  Liquid  saturation  distribution  at  cathode  GDL  for 
isothermal  model  at  Vcen  =  0.55  V  for  AFCL  case. 


A  two-dimensional  simulation  is  performed  when  both  anode  and 
cathode  inlet  gas  streams  are  fully  humidified  ( AFCF)  and  when  the 
anode  inlet  flow  is  fully  humidified  with  low  cathode  flow  humidi¬ 
fication  (AFCL).  An  inlet  stoichiometric  ratio  of  2  is  chosen  for  both 
anode  and  cathode  sides  based  on  a  reference  current  density  of 
1 A  cm-2,  therefore  the  flow  rates  are  fixed  for  anode  and  cathode. 
The  temperature  and  pressure  are  80  °C  and  1.5  atm  at  the  inlets 
of  both  anode  and  cathode  and  the  physicochemical  and  transport 
properties  are  listed  in  Table  4. 

The  numerically  predicted  liquid  saturation  distribution  at  GDL 
cathode  side  at  0.55  V  is  shown  in  Fig.  4a  and  b  for  AFCL  case,  for 
both  the  isothermal  and  non-isothermal  models.  At  this  voltage, 
the  water  may  condense,  but  this  does  not  occur  until  the  water 
vapor  concentration  in  the  gas  reaches  the  saturation  level.  The 
condensation  with  its  location  directly  related  to  the  local  tem¬ 
perature  and  condensation  front  is  then  pushed  downstream.  In 
the  non-isothermal  model,  up  to  a  certain  point  along  its  lengths 
(y  =  0.045  m),  the  GDL  is  free  of  any  liquid  water  after  which  the 
stream  of  liquid  water  starts  to  appear.  In  the  isothermal  model, 
the  water  condensation  and  accumulation  occur  further  upstream 
towards  the  inlet  of  the  channel  at  y  =  0.032  m,  where  the  amount 
of  liquid  saturation  is  higher  in  the  isothermal  compared  with  the 
non-isothermal  model.  This  is  because  in  non-isothermal  model, 
the  temperature  within  the  distance  between  y  =  0.032  and  y  =  L  is 
about  1.8  °C  higher  than  that  of  the  isothermal  model.  The  increase 
in  temperature  causes  an  increase  in  saturation  concentration,  so 
less  water  vapor  changes  phase. 

The  liquid  saturation  level  for  isothermal  and  non-isothermal 
models  for  the  AFCF  case  at  cathode  GDL/catalyst  layer  interface  is 
shown  in  Fig.  5.  In  both  cases,  the  maximum  liquid  saturation  occurs 
near  the  inlet  and  decreases  in  the  flow  direction  due  to  decreasing 
water  production  in  the  channel  length.  The  amount  of  liquid  satu¬ 
ration  along  the  channel  in  the  isothermal  model  is  more  than  the 
non-isothermal  one,  where  the  maximum  deviation  for  both  cases 
appears  at  the  entrance  where  the  cell  temperature  is  highest.  The 
maximum  level  of  liquid  saturation  is  12%  and  10%  for  isothermal 
and  non-isothermal  models,  respectively.  Further  along  the  chan¬ 
nel  the  temperature  drop  causes  the  deviation  to  decrease  as  well 
where  the  vapor  saturation  concentration  is  a  strong  function  of 
temperature. 

In  a  non-isothermal  two-phase  zone  is  water  transport  in  vapor 
diffusion  mode  due  to  variation  in  saturation  vapor  concentra¬ 
tion  with  temperature.  This  mode  of  vapor  diffusion  is  from  the 
high  temperature  to  low  temperature  regions.  For  investigated  of 
this  mode  of  vapor-phase  diffusion  as  driven  by  the  temperature 


Fig.  5.  Liquid  saturation  profiles  at  cathode  GDL/catalyst  layer  interface  along  the 
channel  direction  for  non-  and  isothermal  models  at  Vcen  =  0.55  V  for  AFCF  case. 


gradient,  the  liquid  saturation  level  for  the  AFCF  case  at  cathode 
GDL/catalyst  layer  interface  with  and  without  considering  effect  of 
temperature  gradient  on  vapor-phase  diffusion  is  compared  and 
presented  in  Fig.  6.  When  we  consider  effect  of  temperature  gra¬ 
dient  on  vapor-phase  diffusion,  water  transport  from  the  higher 
temperature  to  the  lower  temperature  regions  (see  Fig.  7).  In  this 
state,  liquid  saturation  slightly  higher,  because  of  the  vapor-phase 
diffusion  that  moves  water  vapor  from  the  water  vapor  production 
in  catalyst  layer  and  the  water  that  income  from  cathode  channel 
to  the  cathode  GDL  for  condensation. 

3.3.  Effect  of  phase  change  on  cell  temperature 

The  temperature  contours  for  a  PEM  fuel  cell  at  0.55  V  for  AFCL 
case  is  shown  in  Fig.  7,  where  the  temperature  rise  is  due  to  the  heat 
generation  caused  by  the  exothermic  nature  of  reactions,  entropic 
heat,  Joule  heating  and  the  water  phase  change.  The  maximum 


Fig.  6.  Liquid  saturation  profiles  at  cathode  GDL/catalyst  layer  interface  along  the 
channel  direction,  with  (state  A)  and  without  (state  B)  considering  vapor-phase 
diffusion  by  the  temperature  gradient,  Vcen  =  0.55  V  for  AFCF  case. 
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Fig.  7.  Temperature  distribution  at  Vcen  =  0.55  V  for  AFCL  case. 


temperature  occurs  near  the  entrance  because  most  of  the  heat 
is  generated  due  to  electrochemical  reaction  at  the  catalyst  layer 
where  the  local  current  density  is  the  highest  at  the  inlet.  The  max¬ 
imum  temperature  difference  is  observed  to  be  3.7  °C  near  the  inlet 
and  at  the  interface  between  membrane  and  catalyst  layer  in  the 
cathode  side  where  major  heat  generation  takes  place.  Unlike  chan¬ 
nels,  velocity  in  the  membrane-electrode  assembly  is  very  low, 
therefore  heat  generation  in  the  cell  is  carried  downstream  smaller 
by  bulk  motion  and  is  primarily  removed  through  the  gas  diffusion 
layer  by  conduction.  This  process  could  be  controlled  by  the  GDL 
thermal  conductivity  or  feed  gas  relative  humidity,  as  inlet  humidi¬ 
fication  strongly  affects  the  degree  of  overall  membrane  hydration. 

The  vapor-phase  diffusion  in  the  non-isothermal  two-phase 
zone  could  cause  to  phase  change  heat  transfer  by  water  evapo¬ 
ration  at  the  hotter  regions  and  vapor  condensation  on  the  cooler 
regions  of  the  cell.  This  mode  of  phase  change  heat  transfer  is 
conventionally  referred  to  as  the  heat  pipe  effect.  Fig.  8  shows 
comparison  the  temperature  profiles  along  the  channel  direction 
at  cathode  GDL/catalyst  layer  interface,  with  and  without  consider¬ 
ing  phase  change  heat  transfer  due  to  vapor-phase  diffusion.  It  can 
be  seen  that  the  first  state  has  a  slightly  smaller  temperature  rise  in 
the  cell.  The  maximum  temperature  deviation  between  two  states 
is  about  0.1  K  that  appears  near  the  inlet. 

To  have  a  better  understanding  of  the  effects  phase  change  on 
temperature  distribution,  a  detailed  comparison  between  the  single 
and  two-phase  models  for  AFCL  and  AFCF  cases  has  been  made. 
The  temperature  profiles  across  the  cell  are  shown  in  Fig.  9  at  three 


Fig.  9.  Comparison  temperatures  profiles  across  cell  direction  at  inlet,  mid  and  outlet 
of  channel  at  Vcen  =  0.55  V  between  single-  and  two-phase  models  for  AFCL  case. 


different  sections,  these  are  the  inlet,  middle  and  outlet  zones  for 
the  AFCL  case  for  the  single  and  two-phase  models.  In  the  two- 
phase  model,  the  phase  change  dose  not  occur  at  the  beginning  of 
the  channel  and  the  heat  release  due  to  phase  change  would  have  a 
smaller  effect  on  the  maximum  temperature  that  appear  in  the  cell 
inlet.  The  maximum  difference  occurs  near  the  exits  of  both  the 
cathode  channel  and  the  GDL  for  both  the  single  and  two-phase 
models,  where  phase  change  takes  place  and  the  heat  generated  by 
condensation  would  increase  the  temperature. 

For  the  AFCF  case,  the  cell  temperature  at  all  zones  for  two-phase 
model  is  higher  due  to  the  phase  change  that  appear  at  the  inlet 
as  shown  in  Fig.  10.  The  maximum  deviation  in  temperature  for 
both  models  occurs  near  inlet  at  GDL  cathode  and  cathode  channel. 
The  temperature  distribution  is  seen  to  be  in  good  agreement  with 
the  available  data  [32]  where  the  minor  difference  could  be  due  to 
single-phase  assumption  versus  our  two-phase  model. 


Fig.  8.  Temperatures  profile  along  the  channel  direction  at  cathode  GDL/catalyst 
layer  interface,  with  (state  A)  and  without  (state  B)  considering  phase  change  heat 
transfer  by  vapor-phase  diffusion,  Vcen  =  0.55  V  for  AFCF  case. 


Fig.  10.  Comparison  temperatures  profiles  across  cell  direction  at  inlet,  mid  and 
outlet  of  channel  at  Vcen  =  0.55  V  between  single-  and  two-phase  models  for  AFCF 
case. 
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Fig.  11.  Comparison  of  polarization  curves  from  isothermal  and  non-isothermal 
calculations. 


3.4.  Comparison  of  isothermal  and  non-isothermal  model 
predictions 

A  two-phase  non-isothermal  transport  phenomena  in  PEM  fuel 
cell  operation  is  studied  in  detail  and  comparison  between  the 
isothermal  and  non-isothermal  model  predictions  for  the  AFCL 
case,  has  been  made.  A  complete  steady  state  polarization  curve 
predicted  by  the  isothermal  and  non-isothermal  models  are  shown 
in  Fig.  11.  At  low  current  densities,  the  cell  voltage  prediction  for 
the  non-isothermal  model  is  higher  than  that  of  the  isothermal 
model.  This  is  because  in  non-isothermal  model,  the  cell  temper¬ 
ature  would  increase  the  exchange  current  density  as  shown  in 
Table  3.  At  higher  current  densities  (0.3  </ave  <  1.2  A  cm-2),  the  cell 
voltage  would  be  higher  for  the  non-isothermal  model,  due  to  the 
fact  that  the  increase  in  cell  temperature  would  improve  the  cell 
potential  owing  to  the  reduced  losses  in  the  cell  and  increase  in 
the  ionic  conductivity  that  leads  to  smaller  resistive  loss  in  the 
membrane.  This  can  also  reduce  transport  and  activation  losses. 
When  the  average  current  density  is  more  than  1.2  A  cm-2,  the  cell 
operates  at  a  limited  mass  transfer  regime  and  there  is  a  signifi¬ 
cant  amount  of  water  generated  due  to  electrochemical  reaction 
in  the  cathode,  leading  to  liquid  water  formation  and  flooding  of 
the  cathode  GDL.  In  the  isothermal  model,  saturation  concentra¬ 
tion  is  smaller  due  to  decreasing  temperature,  resulting  in  higher 
amounts  of  liquid  water  formation  in  the  cathode  GDL.  Therefore, 
the  cell  performance  is  hampered  in  comparison  with  the  non- 
isothermal  model.  However,  high  current  densities  could  increase 
water  vapor  partial  pressure,  enhancing  mass  transport  related 
losses. 

4.  Conclusions 

A  two-dimensional,  two-phase,  non-isothermal  model  where 
the  conservation  equations  are  coupled  with  electrochemical 
kinetics  has  been  developed  for  a  PEM  fuel  cell,  based  on  a  multi¬ 
phase  mixture  formulation  to  investigate  the  interaction  between 
heat  and  water  transport  together  with  phase  change.  The  main 
focus  was  to  assess  the  results  from  coupling  heat  and  trans¬ 
port  phenomena  in  the  fuel  cell.  The  following  conclusions  are 
drawn: 


•  The  temperature  distribution  is  the  most  important  parame¬ 
ter  affecting  the  two-phase  water  transport,  which  controls  the 
amount  and  location  of  liquid  saturation. 

•  Partial  humidification  of  the  inlet  air  causes  the  liquid  water  to 
build  up  down  stream  at  cathode  GDL,  whose  generated  location 
is  temperature  sensitive,  so  in  the  non-isothermal  model,  where 
temperature  increase  enhances  saturation  concentration,  a  larger 
portion  along  the  channel  length  will  be  free  of  liquid  water. 

•  Full  humidification  of  the  air  at  the  inlet  results  in  condensation 
and  flooding  at  full  face  of  GDL  cathode,  where  the  maximum  liq¬ 
uid  saturation  at  GDL  cathode  moves  further  up  stream  towards 
the  inlet  and  decreases  in  the  flow  direction  along  the  chan¬ 
nel.  The  amount  of  liquid  saturation  in  the  isothermal  model  is 
more  than  the  non-isothermal  model,  and  the  maximum  devia¬ 
tion  occurs  at  the  entrance  due  to  the  extremely  high  temperature 
at  the  inlet. 

•  The  PEM  fuel  cell  temperature  distribution  is  highly  voltage  sensi¬ 
tive;  this  model  predicts  the  temperature  differential  to  be  about 
3.7  °C  at  a  chosen  voltage  of  0.55  V.  The  maximum  temperature 
occurs  near  the  inlet  at  the  interface  between  the  membrane  and 
catalyst  layer  on  the  cathode  side  where  major  heat  generation 
takes  place. 

•  At  partial  humidification  of  inlet  air,  phase  change  would  have  a 
smaller  effect  on  the  maximum  temperature  that  appears  in  the 
cell  inlet,  also  the  maximum  temperature  difference  for  single 
and  two-phase  models  occur  further  down  stream  towards  the 
exit  of  cathode  channel  and  its  GDL. 

•  At  full  humidification  of  the  inlet  air,  the  cell  temperature  at  all 
regions  for  the  two  phase  model  is  higher  due  to  the  phase  change 
that  occurs  at  the  inlet,  where  the  maximum  temperature  devia¬ 
tion  for  both  models  occurs  further  up  stream  in  cathode  channel 
and  its  GDL. 

•  Vapor-phase  diffusion  by  providing  a  new  mechanism  for  heat 
removal  from  the  cell  affects  the  temperature  distribution. 
Instead,  water  transport  via  vapor-phase  diffusion  due  to  the  tem¬ 
perature  gradient. 

•  The  accuracy  of  performance  prediction  by  the  non-isothermal 
model  is  better  than  the  isothermal  one,  in  which  the  temper¬ 
ature  rise  may  increase  the  exchange  current  density  and  ionic 
conductivity  or  reduce  flooding  in  the  fuel  cell. 

•  The  present  study  successfully  demonstrates  the  importance  and 
accuracy  of  a  coupled,  two-phase  heat  and  water  transport  model 
in  a  PEM  fuel  cell. 
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